On the Convergence of Collocation Methods for Boundary Integral Equations on Polygons

The integral equations encountered in boundary element methods are frequently solved numerically using collocation with spline trial functions. Convergence proofs and error estimates for these approximation methods have been only available in the following cases: Fredholm integral equations of the second kind [4], [7], one-dimensional pseudodifferential equations and singular integral equations with piecewise smooth coefficients on smooth curves [2], [3], [17], [26)-[29], and some special results on the classical Neumann integral equation of potential theory for polygonal plane domains [5], [8], [9]. Here we give convergence proofs for collocation with piecewise linear trial functions for Neumann's integral equation and Symm's integral equation on plane curves with corners. We derive asymptotic error estimates in Sobolev norms and analyze the effect of graded meshes. 0. Introduction. In this paper we give convergence proofs and asymptotic error estimates in Sobolev norms for collocation with piecewise linear spline trial functions applied to two basic integral equations of potential theory on plane polygons, namely the integral equation of the second kind with the double layer potential (" Neumann's integral equation"), and the integral equation of the first kind with the simple layer potential ("Symm's integral equation"). We use an idea of Arnold and Wendland [2], namely considering Dirac delta functions (the " test functions" in the collocation method) as second derivatives of piecewise linear functions. Therefore, similar results as presented here should be possible for splines of higher odd order. Corresponding results for even-order splines are not yet available. Thus, for one of the simplest methods of numerically solving Dirichlet's problem on a plane domain with corners, the midpoint collocation with piecewise constant trial functions for the first-kind integral equation with the simple layer potential, convergence is still an open problem. The method of Fourier series that yields the convergence proof in the case of a smooth boundary [27] cannot be applied in the presence of corners. We apply the method of local Mellin transformation that has previously been used to derive error estimates for Galerkin methods for a wide class of operators, including those occurring in boundary element methods in acoustics, electromagnetism, and elastostatics [U]-[14]. Thus, it is to be expected that also the techniques presented here will apply to a rather large class of integral equations. For example, Received June 9, 1986. 1980 Mathematics Subject Classification (1985 Revision). Primary 65N35, 65R20; Secondary 65D07, 45L10. *This author was supported by the National Science Foundation under Grant DMS-8501797 and Grant DMS-8603954. 461 ©1987 American Mathematical Society 0025-5718/87 $1.00+ $.25 per page License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use 462 MARTIN COSTABEL AND ERNST P. STEPHAN the case of singular integral equations of Cauchy type will be treated in a forthcoming paper. Let T be a connected closed plane curve composed of smooth arcs TJ, j = 1,..., J, that meet'at the corner points Zj at interior angles Uj e (0,27r). The Sobolev spaces HS(T) are defined for s > 0 being the restriction of HS+1/2(R2) to I\ for s < 0 by duality: HS(T) = H-S(T)', and H°(T) = L2(T). It is known [11], [12] that for |s| < 3/2, the space HS(T) may equivalently be defined as the corresponding Sobolev space on the arc length parameter interval, transferred to T by the parameter representation map. We consider the following two integral equations on T: (0.1) (l + K)u=f, (0.2) Vu = f. Here the operator K of the double layer potential is defined by (0.3) Ku(z):= -];f u(S)-r^-rlog\z-!;\ds{n, ir Jr on(S) where s(f ) is the arc length on T and 3/3«(f ) denotes the derivative with respect to the normal vector at f e T pointing into the interior of T. The operator V of the simple layer potential is defined by (0.4) Fii(f)--( MU)log|z-n aj G (i> 1) by v= min{~* 2ttu}' a°:= mm{aj\J = h---,J}Similarly, V: HS(T) -» HS+1(T) is continuous and bijective for all í e (\ a0, — 2 + ao)> provided the analytic capacity of T is not equal to one. We shall assume this in the sequel. For the collocation method, we need a grid A^ = {xx,...,xN} c T, the x, being both the collocation points and the meshpoints of the trial functions. By Sl(AN) we denote the N-dimensional space of splines of order 1, i.e., each u e S1(hN) is a continuous function on T that is a linear function of the arc length on each of the segments xnx v n = 0,..., N 1, where x0:= xN. Let h:= max{|jc„ + 1 x„||n = 0,...,N 1}. We do not impose a uniformity condition on AN, but assume only that h -* 0 as TV tends to infinity. For the second-kind integral equation (0.1), the collocation method is the following: Find uN e Sl(AN) such that (0.5) (l + K)uN(xn)=f(x„), n = l,...,N. For the first-kind integral equation (0.2), the collocation equations are (0.6) VuJx„) = f(xn), n = l,...,N, License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use COLLOCATION METHODS FOR BOUNDARY INTEGRAL EQUATIONS 463 but we shall have to modify (0.6) in some way in order to obtain a convergence proof (see (3.3) and (3.6)). The paper is organized as follows: In Section 1 we present some facts on the convergence of general projection methods. They are stated in a form which is convenient for the application to collocation methods and allow easy incorporation of compact perturbations as well as localization arguments. In Section 2 we prove convergence and stability in the Hl Sobolev norms for the approximation scheme (0.5) for the second-kind integral equation (0.1). In Section 3 convergence and stability results for two modifications of the scheme (0.6) for the first-kind integral equation (0.2) are shown. In the final Section 4 we investigate the asymptotic orders of convergence. For the case of the first-kind integral equation, where we have to use weighted Sobolev norms, we prove a new approximation result and we show that the use of suitably graded meshes yields convergence of the same order as for smooth curves. 1. On the Convergence of Projection Methods. We need some results on the convergence of projection methods, including compact perturbations and spaces with two norms. Such results are well known [16], [19], [24], but we present a formulation that is particularly adapted to the present case. As the lemma in question might be of independent interest, we also include a complete proof. Let X and Y be Banach spaces and A : X -> y be bijective and continuous. For the approximate solution of the equation (1.1) Au=f we assume that we have a sequence of finite-dimensional subspaces VN czX, TNc r, dimJ^ = dimTN < oo (N e N) and we replace Eq. (1.1) by the relation for uN G VN, (1.2) (t,AuN) = (t,f) for allí G 7„. Here the brackets denote the duality between the space Y and its dual Y'. We make the following assumptions: (i) There exist bounded operators PN: Y' -> TN that converge on Y' strongly to the identity operator. (ii) There is a Banach space X0, continuously embedded in X (hence, ||x||^< C||jc|| % for all x g X0 and some constant C). (iii) For all N there holds VN c X0. (iv) For all N we are given a mapping QN: VN -* TN and a constant M such that (1.3) \{QNv,Aw)\^M\\v\\x\\w\\Xu for all u g VN, we X0, /V g N. (v) There exists a collectively compact sequence of operators CN: X -> X' in the sense of [1] and a constant y > 0 such that (1.4) \(QNv,Av) + (CNv,v)\>y\\v\\2x for all o g VN, N g N. Lemma 1.1. Under the above conditions (i)-(v) there exists N0 g N such that for all N ^ N0 the system (1.2) has a unique solution uN g Vn for any f g Y. There is a constant C such that for this "approximate solution" uN and the "true solution" u License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use 464 MARTIN COSTABEL AND ERNST P. STEPHAN there holds (1.5) IKIIa- N0, (1-6) II« uN\\x < Cinf{||« ö||x0|ö g VN). Proof. For abbreviation, we write || • ||:= || • H^ and || • ||0:= || • H^. The derivation of the quasi optimality (1.6) from the stability (1.5) and the unique solvability of the system (1.2) is standard: Denote the solution operator u »-» uN by GN. Then GN: X0 -» (VN, || • ||) is a projection operator. Its norm is bounded by C for all N > 7V0 by (1.5). Thus for all VN: \u ~ un\\ =||« ö GN(u «)|| ^ ||u ü|| + C||w w||0. Hence the assertion (1.6) follows. For the proof of unique solvability of (1.2) and stability estimate (1.5) we consider first the special case where all operators CN in assumption (v) vanish. (Actually, this is not so special: The existence of QN with (1.4) for CN = 0 is also necessary for stability, cf. [17].) Then from (1.4) it follows that the solution of (1.2) is unique, namely: If (t, Av) = 0 for all t G TN and some v g Vn, then yII!;II <<\{QnV,Au)\0, hence v = 0. As (1.2) is represented, after choosing bases in VN and TN, by an N X N system of equations, the existence of uN follows from uniqueness. The stability estimate (1.5) follows from (1.2), (1.3), (1.4): Il M < -\(QNuN,AuN)\ = -\{QNuN,Au)\ < —1| uN || || u ||0. Now we consider the general case with nonvanishing perturbations CN. We define new operators QN and show that all assumptions are satisfied for large /V if we replace QN by QN and CN by 0. Thus we reduce the general case to the special case considered above. Define QN:= QN + PNA' lCN =QN + A'~lCN-{\ PN)A"lCN. Here A'_1 is the inverse of the isomorphism A': Y' -» X' adjoint to A. The norms of PNA'~lCN: X ^> T'are bounded, hence \(QNv,Aw)\<\(QNv,Aw)\ +\(pNA'-lCNv,Aw)\ < M\\v\\ \\w\\0 + Mx\\v\\ \\w\\ « (M + CMX)||v|| ||w||0. Thus (1.3) holds for QN and all N. As 1 PN -> 0 strongly on Y' and the operators A'~lCN: X -+ Y' are collectively compact, A'(l PN)A'~lCN: X -> X' tends to zero in operator norm. If we denote the operator norm of A'(l PN)A'~lCN by 8N, License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use COLLOCATION METHODS FOR BOUNDARY INTEGRAL EQUATIONS 465 we obtain \(QNv,Av)\ =\(QNv,Av) + (A'-lCNv,Av) + {{I PN)A'-lCNv, Av)\ = \(QNv,Av) + (CNv,v) + (A'(l PN)A'-lCNv,ü)\ >(y-8N)\\v\\2. Thus, if N is large enough to imply ô^ < y, the corresponding estimate (1.4) holds for QN, and CN replaced by zero. This completes the proof. D Remark 1.2. We shall need the lemma only for the case of QN and CN not depending on /V. Thus QN = Q: X0 -> Y' will be a linear operator satisfying (1.7) QVn^Tn for all TV g N, and CN = C: X -* X' will be a compact operator, or equivalently, the quadratic form v -* (Cv, v) appearing in (1.4) will be completely continuous on X. Remark 1.3. The operators PN: Y' -» TN are not explicitly needed for (1.2). Only their existence is used in the proof. If Y' is a Hilbert space, we can take the orthogonal projections onto TN. We then must assume that TN -* Y' in the sense that for all / G Y' there is a sequence tN G TN converging to t. By duality and the reflexivity of Y, this can be formulated as the following condition. , . If ve Fand lim N^x{tN,y) = 0 for each sequence tN g Tn, (1'8) thenv-0. We shall use this condition later on instead of (i) above. The Gârding type inequality (1.4) can be localized by means of a partition of unity. We formulate this result for the situation of spaces with two norms but with QN and CN not depending on N. Thus, we make the following assumptions: Q: X0 -» T' is a linear operator with A'Q: X -» X¿ bounded (according to (1.3)). There exist bounded linear commuting operators a(_/' = 1, ...,m) on X and b¡ (j = 1,..., m) on Y such that (a)Zfmlaj 1 oa X; (ß) Bf:= AOj bjA is compact from X to AX0 c Y; i.e., A~lBf: X -> X0 is compact; (Y) B} := Qoj b'jQ is compact from X to Y'\ (8) For every k = l,..., m there is a compact operator Ck: X -> X' and a constant y^ > 0 such that (1.9) Re((Qakv,Aakv) + (Ckv,v)) > yk\\akv\fx for all o G X. Lemma 1.4. Let the assumptions (a)-(8) be satisfied. Then there exists a compact operator C: X -» X' and a constant y > 0 such that (1.10) Re((Qv,Av) + (Cv,v))>y\\v\\2x forallv&X.; In particular, (1.4) holds. Proof. From (ß) and (y) it follows that for j, k = 1,..., m, (Qa2v,Aa2v) = (QakajV, AakajV) + (Ckjv,v) License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use 466 MARTIN COSTABEL AND ERNST P. STEPHAN with Ckj: X -* X' compact. Here one has to use that Bj"Q: X -> X' is compact. There follows (Qu,Av) = E \Qa2v,Aa2kv j,k = \ = E {{QakajV,Aakajv) + (Ckjv,v)). j,k=i By (1.9), Re{(QakajO,Aakajv) + (CkajV,ajv)) > yk\\akajv\\ , hence m Re(Qv,Av)> E (yJIûa^II Re(ckjv,v) Re(ckajV,ajV 7,* = 1 2 > -yl| t71| Re(Ct;,i;) with y := ¿min{ yj /c = 1,..., m} and C = EM(Ct> + a'£ka¡). D Remark 1.5. Note that in this formulation the finite-dimensional spaces K^ and TN do not appear. Therefore, this local principle is very easy to apply. Compare also the local principles of Prössdorf [22] and Silbermann [17]. 2. The Second-Kind Integral Equation. In order to show convergence of the collocation scheme (0.5) for the integral equation (0.1), we apply Lemma 1.1 to the following situation: X0 = X = Y = Hl(T); A = 1 + K; VN:= S\àN); TN:= S-\LN),vhsn (2.1) S-l{AN):=spm{8(x-xn)\n = l,...,N}. We have to check the assumptions of Lemma 1.1. The abstract Galerkin equations (1.2) coincide with the collocation equations (0.5). Furthermore, assumption (i) of Lemma 1.1 is satisfied in view of Remark 1.3. Namely, Hl(T) is continuously embedded into C(T) and the condition h -> 0 implies that every point f g T is an accumulation point of a sequence ^ e AÄ. Thus the hypothesis in condition (1.8) is only satisfied if y = 0 on I\ Next we define (2.2) Q:=D2, i.e., the second distributional derivative with respect to the arc length. Then clearly (2.3) QS^A^cS-^A») holds, i.e., (1.7) is satisfied. Note that our lemma does not require QN to be bijective! (The latter property is frequently assumed in other approaches [2].) We have to show the two estimates (1.3), (1.4) in the case QN = Q. The first one follows by continuity: Let v, w G H1^); then |<ßi;,(l + K)w)\ =\-(Dv,D(l + K)w)LHV)\ 0 anda compact operator C: Hl(T) -* //_1(F) such that (2.5) \(Dv,D(l + K)v) + (Cv,v)\>y\\v\\2HHr) for all v g Hl(T). For the proof of (2.5) we use a partition of unity and Lemma 1.4 to reduce (2.5) to the corresponding estimate on a reference angle. Let Tw = e""R + U R+ be this reference angle. If we use Tu to parametrize a neighborhood of one of the corners z-, then the operator induced by K differs from the operator of the double layer potential defined on Ta only by an operator that is compact on Hl; see [10]. Thus we only need to consider the case that T and Tu coincide on a neighborhood of the origin, and K is defined on Tu. We then have to show Lemma 2.2. There is y > 0 such that (2.6) Re(Dv,D(l + K)v) > y\\vÙiTJ for all v G //^r^) with support in a fixed compact set. The constant y may depend on this compact set and the angle u, but not on v. Proof. We proceed analogously to the proof of Gàrding's inequality for 1 + K in [11]. The operator K maps even and odd functions on Tu to even and odd functions, respectively. Therefore, it suffices to show (2.7) Re(Dv,D{l±Kjv)L2(Rt)>y\\v\\2HHR+) for all v G C0°°[0, oo) with Kj(x):Iflra(—J-tau')«/»'. By the Parseval relation for the Mellin transform we obtain Re(Dv,D(l ± Kjv) (2.8) = ^-Re f \X + i\2 2ir hm\=-\/2 I sinh(7T -to)(\ + /')\l.,. ., .2 ,, \ sinh7r(\ + /) /' ' Here the Mellin transform is defined by J/-00 x'x-lv(x)dx, 0 and we use Dv(X)= -i(X + i)v(X + i), K^(X)=-i$inh{\-")Xv(X) forImAG(-l,l)[ll]. " sinhwA License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use 468 MARTIN COSTABEL AND ERNST P. STEPHAN Now we have sinh(7T u)X sinhirX ir sin — = :1 — q < 1 for all À with Im À = — (In [11] we used this estimate for Im A = j.) Therefore, we can estimate (2.8) from below by 1 — qC ,-, .|2,A/, .\i2 ,, . ... „2 „ „2 f \X + i\2\v(X + i)\ dX=(l-q)\\Dv\\2L2iR+)>y\\v\\2Hi(R+). lir J\m\=-\/2 D Proof of Lemma 2.1. Choose a partition of unity x, e CTM(R2) with ET-iX/ = 1 on T such that Xj• > 0 and x/= vX/ e 0° and such that the support of Xj contains exactly one corner Zj. (Thus m = J.) Then a3 = è= ¿»j in Lemma 1.4 is the operator of multiplication by x¡Then the commutator of A with a} is compact on HX(Y) (actually it maps L2(T) into H1^) continuously), and the commutator of Q with üj maps Hl(Y) into L2(T) and hence is compact from HX(Y) to H~l(T). Together with Lemma 2.2, we see that all assumptions of Lemma 1.4 are satisfied, and its conclusion (1.10) gives (2.5). □ We can now apply Lemma 1.1 and find immediately Theorem 2.3. There is an N0 g N such that for all N ^ N0 and all f G H\T) the system of collocation equations (0.5) has a unique solution uN G S^A^,). There is a constant C independent of N and u such that ll"jvll//'(r)< C||i/||Hi(r) and II"«Jlff'tn^ Cinf{||nù||Hi(r)|ùG S^A/v)}. Remark 2.4. Note that the grids AN need not be uniform. Nor do they have to include the corner points. One can conjecture that the estimates in Theorem 2.3 hold for arbitrary bounded Lipschitz domains, probably also in higher dimensions. Up to now, however, there exists no proof avoiding the use of Mellin transformation. 3. The First-Kind Integral Equation. We consider the collocation scheme (0.6) for the integral equation (0.2). The natural choice Q = D2; X = H+l/2(T), A = V is not useful, because in general then (Qv, Av) = oo, i.e., Q does not map into y = (AX)' c H~3/2(Y). We present two variations of this natural choice. Let us first define H^2{Y):= {«G//1/2(r)|ù/G//1/2(r)(;' = l,...,/), (3.1) where üj — u on TJ, ùJ■ = 0 on T \ TJ j. The norm in Hl/2(T) is II2 •— V II II2 llH1/2(r)-Zw \\uj\\nl/2(ry 7 = 1 Then Hl/2(T) is the completion of C*(T\{zx,..., Zj}) in this norm. Hl/2(T) is densely embedded in H1/2(T). The operator V maps Hi/2(T) continuously into H3/2(T) but it is not surjective, as the //"1/2(r)-solution u of Eq. (0.2) with smooth License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use COLLOCATION METHODS FOR BOUNDARY INTEGRAL EQUATIONS 469 / in general behaves like 0(\z z¡\tt'~x) near z-, where ay was defined in the introduction [11]. Thus, in general, u G HS(T) for s > a0 \ G (0, \). We define (3.2) S(AN):= {v£Sl(AN)\v(Zj) = 0,j = l,...,j}, and we assume that {zx,...,zj} cAN. Thus, dim S(AN) = N J. Now choose J functions ijx,..., tj7 g H3/2(T) such that r\j(zk) = 8jk (j,k = l,...,J). Define the projection operator R: H3/2(T) -» //3/2(r)n^1/2(r)by tfg(z):=g(z)I g(z,.)ij,(z). 7 = 1 Then the adjoint operator R' acts in S_1(^/v) as follows: j R'8(z xk) = 8(z -xk)E 1j(xk)8(z zj) (xk g An). 7 = 1 If we take the (N J )-dimensional range space of R' as space TN of test functions, we arrive at the modified collocation equations j j VuN{xn)E VuN(zj)t\j(x„) =f(xn)¿Zf(zj)flj(x„) (3.3) 7 = 1 7 = 1 (n = l,...,N),uN£S(AN). These are satisfied by solutions in S(AN) of (0.6) but not conversely. Now we have the following situation: XQ = X= H1/2(T); A = V; Y = AX; VN = S(AN); TN = R'Sl(AN); Q = R'D2. Furthermore, we can assume that tj G Y for j' = l,...,J. This implies (1.8). Then it is easily seen that the assumptions (i)-(iv) of Lemma 1.1 are satisfied. It remains to show the inequality (1.4). Lemma 3.1. There is a compact operator C: Hl/2(T) -* //1/2(r)' and a constant y > 0 such that Re((Dv,DRVv) + (Cv,v)) > y\\v\\\i/HT) for all v g Hl/2(T). Proof. By localization to the reference angle Ta, decomposition into even and odd parts, and density arguments, we see that we have to show (3.4) Re{(Dv,D(V0± Vjv) + (CM» > yM|W+> for all v g Cq°(R+) with support in a fixed compact set. Here, (3.5) VMx)-= --/ logl -e""" Hy)dyir J0 y License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use 470 MARTIN COSTABEL AND ERNST P. STEPHAN Mellin transformation gives [11] ?3W-2fSa**<*-'> 0"* «<•.»>• By the Parseval relation we obtain (Dv,D(V0± Vjv) = xf \X + i\2i>(X + i) ¿IT /uu _i n |ß(X)</A. 'lmX= -1/2 coshw(A + /') + cosh(flu)(X + i) (X + i) sinhw(A + i) For Im À = \ we have v(X + i) = v(X) and X + i = X. Shifting the path of integration to Im X = 0, we thus obtain (Dv,D(V0±Vjv) 1 r C0Sh77(A + l') ± cosh(77 u)(X + i) , . , . ,2 = "t— / A-~,-7^-TTÜ(A) «A 2"'-/ImX-0 Sinh77(A + /) -j-j m+(X)\v(X)\ dX,